/*

cimahis - Software para el analisis de imagenes de corte histologico

Copyright (C) Andrés De Abreu 2010 <jadabreu @ gmail . com>


Este programa es software libre. Puede redistribuirlo y/o modificarlo bajo los términos
de la Licencia Pública General de GNU según es publicada por la Free Software Foundation,
bien de la versión 2 de dicha Licencia o bien (según su elección) de cualquier versión
posterior.

Este programa se distribuye con la esperanza de que sea útil, pero SIN NINGUNA GARANTÍA,
incluso sin la garantía MERCANTIL implícita o sin garantizar la CONVENIENCIA PARA UN
PROPÓSITO PARTICULAR. Véase la Licencia Pública General de GNU para más detalles.

Debería haber recibido una copia de la Licencia Pública General junto con este programa.
Si no ha sido así, escriba a la Free Software Foundation, Inc.,
en 675 Mass Ave, Cambridge, MA 02139, EEUU.

    //ACERCA DE histograma.cpp
    *
    * Implementacion de la clase histograma y de las subclases histogramaTinte
    * e histogramaEstandar.
    *

*/


#include <iostream>
#include <opencv/cv.h>
#include <opencv/highgui.h>
#include "color/Histograma.h"
#include "math.h"
#include "vector"
#include "stdio.h"

using namespace std;


///////////////////////Constructores de las clases concretas///////////////////
histogramaTinte::histogramaTinte()
{
    hist=0;
    nBarras=360;
    media=0;varianza=0;
    
    float rangosCanal[] = {0,180};  // internamente OpenCV representa los valores de tinte de 0 a 180 debido a que en 8 bits solo caben 255 valores, optarnon por representar el tinte como 360/2
    float* rangos[] = {rangosCanal};
    
    hist=cvCreateHist(1, &nBarras, CV_HIST_ARRAY, rangos, 1); 
  

}

histogramaEstandar::histogramaEstandar()
{
    hist=0;
    nBarras=256;
    media=0;varianza=0;
        
    float rangosCanal[] = {0,255}; //en RGB todos los canales tienes rango de valores de 256
    float* rangos[] = {rangosCanal};

    hist=cvCreateHist(1, &nBarras, CV_HIST_ARRAY, rangos, 1); // Funcion de Opencv para crear un histograma por defecto se crea de tipo array, el uno al final es un flag de uniformidad
   
}



/////////////////////////////////Calcular Histograma///////////////////////////////////

void histogramaTinte::calcularHist (IplImage *img)
{
  //funcion de OpenCv para calcular histogramas, se le pasa la imagen, el histograma, un flag para indicar si
    //el histograma es acumulativo y al final la mascara que se desea aplicar..
    cvCalcHist( &img, hist, 0, NULL ); 
   
}

void histogramaEstandar::calcularHist (IplImage *img)
{
    //funcion de OpenCv para calcular histogramas, se le pasa la imagen, el histograma, un flag para indicar si
      //el histograma es acumulativo y al final la mascara que se desea aplicar..
    cvCalcHist( &img, hist, 0, NULL );


}



//////////////////////////////////////////////////////////////////////////////
//Funcion que transforma un valor de tinte a 
//un escalar con los valores rgb correspondientes
//esta funcion es necesaria para graficar los valores del canal H
// en imagenes HSV
CvScalar hsv2rgb( float hue )
{
    int rgb[3], p, sector;
    static const int sector_data[][3]=
        {{0,2,1}, {1,2,0}, {1,0,2}, {2,0,1}, {2,1,0}, {0,1,2}};
    hue *= 0.033333333333333333333333333333333f;
    sector = cvFloor(hue);
    p = cvRound(255*(hue - sector));
    p ^= sector & 1 ? 255 : 0;

    rgb[sector_data[sector][0]] = 255;
    rgb[sector_data[sector][1]] = 0;
    rgb[sector_data[sector][2]] = p;

    return cvScalar(rgb[2], rgb[1], rgb[0],0);
}

///////////////////////////////Dibujar Histogramas//////////////////////////////////

IplImage * histogramaTinte::dibujarHist()
{
	
    IplImage *img;


    img = cvCreateImage(cvSize(360,281),IPL_DEPTH_8U,3); //imagen donde se dibujara el histograma

        int AnchoBarra,ancho,alto;
        float MinVal,MaxVal,valor;

		
	
        cvGetMinMaxHistValue( hist, &MinVal, &MaxVal); //obtener el valor maximo y minimo
        cvConvertScale( hist->bins, hist->bins, MaxVal ? 255. / MaxVal : 0., 0 ); // Ajustar el alto de las barras para que quepan todas en una imagen

	ancho = img->width;
	alto = img->height;
        AnchoBarra= ancho / nBarras;
	
	
    cvRectangle(img, cvPoint(0,0), cvPoint(ancho,alto), CV_RGB(255,255,255),-1);// dibujar el fondo del histograma
	

    for(int i=0; i < nBarras; i++)
    {
            valor = cvRound( cvGetReal1D(hist->bins,i)*alto/255 ); //obtener valor i
            CvScalar color = hsv2rgb(i*180.f/nBarras); //cambiar de HSV a RGB para graficar

            cvRectangle( img, cvPoint(i*AnchoBarra,alto), //dibujar barra
                                 cvPoint((i+1)*AnchoBarra,alto - valor),
                                 color, -2, 8, 0 );
    }
    return img; //devolver la imagen donde se dibujo el histograma

}

IplImage * histogramaEstandar::dibujarHist()
{


    IplImage *img;


    img = cvCreateImage(cvSize(512,281),IPL_DEPTH_8U,3); //imagen donde se dibujara el histograma

    int AnchoBarra,ancho,alto;
    float MinVal,MaxVal,valor,valorAnterior;


    cvGetMinMaxHistValue( hist, &MinVal, &MaxVal); // obtencion de los valores minimos y maximos
   

    cvConvertScale( hist->bins, hist->bins, MaxVal ? 255. / MaxVal : 0., 0 ); // Ajustar el alto de las barras para que quepan todas en la imagen


    ancho = img->width;
    alto = img->height;
    AnchoBarra= ancho / nBarras;


cvRectangle(img, cvPoint(0,0), cvPoint(ancho,alto), CV_RGB(255,255,255),-1); // se dibuja el fondo del histograma

valorAnterior=0;

for(int i=0; i < nBarras; i++)
{
        valor = cvRound( cvGetReal1D(hist->bins,i)*alto/255 ); //obtencion de el valor i
       

        cvRectangle( img, cvPoint(i*AnchoBarra,alto), // dibujamos la barra
                           cvPoint((i+1)*AnchoBarra,alto - valor),
                             cvScalar(0,0,0,0), -1, 4, 0 );
        valorAnterior=valor;
        
}

    return img;

}

CvHistogram * histograma::pasarApuntador()
{
    return hist; //metodo para obtener el apuntador al histograma
}

IplImage * histograma::dibujarGrupo(histograma *h2,histograma *h3) //funcion diseñada para dibujar en la misma imagen 3 histogramas
{
    IplImage *img;
    CvHistogram *hist2,*hist3;

    hist2=h2->pasarApuntador();
    hist3=h3->pasarApuntador();


    img = cvCreateImage(cvSize(768,281),IPL_DEPTH_8U,3);

    int AnchoBarra,ancho,alto;
    float MinVal,MinVal2,MinVal3,MaxVal,MaxVal2,MaxVal3,valor,valor2,valor3,valorAnterior,valorAnterior2,valorAnterior3;


    cvGetMinMaxHistValue( hist, &MinVal, &MaxVal);
    cvGetMinMaxHistValue( hist2, &MinVal2, &MaxVal2);
    cvGetMinMaxHistValue( hist3, &MinVal3, &MaxVal3);


    cvConvertScale( hist->bins, hist->bins, MaxVal ? 255. / MaxVal : 0., 0 );
    cvConvertScale( hist2->bins, hist2->bins, MaxVal2 ? 255. / MaxVal2 : 0., 0 );
    cvConvertScale( hist3->bins, hist3->bins, MaxVal3 ? 255. / MaxVal3 : 0., 0 );



    ancho = img->width;
    alto = img->height;
    AnchoBarra= ancho / nBarras;


cvRectangle(img, cvPoint(0,0), cvPoint(ancho,alto), CV_RGB(255,255,255),-1);

valorAnterior=0;
valorAnterior2=0;
valorAnterior3=0;

for(int i=0; i < nBarras; i++)
{
        valor = cvRound( cvGetReal1D(hist->bins,i)*alto/255 );
        valor2 = cvRound( cvGetReal1D(hist2->bins,i)*alto/255 );
        valor3 = cvRound( cvGetReal1D(hist3->bins,i)*alto/255 );


        cvLine(img,cvPoint(i*AnchoBarra,alto-valorAnterior),cvPoint((i+1)*AnchoBarra,alto-valor),
               cvScalar(255,0,0,0),2,4,0);
        cvLine(img,cvPoint(i*AnchoBarra,alto-valorAnterior2),cvPoint((i+1)*AnchoBarra,alto-valor2),
               cvScalar(0,255,0,0),2,4,0);
        cvLine(img,cvPoint(i*AnchoBarra,alto-valorAnterior3),cvPoint((i+1)*AnchoBarra,alto-valor3),
               cvScalar(0,0,255,0),2,4,0);


        valorAnterior=valor;
        valorAnterior2=valor2;
        valorAnterior3=valor3;

}

    return img;

}



float histograma::obtenerMedia() // calcular la media del histograma
{
	int acumulador;
	float resultado,frecuencia,nElementos;
	nElementos=0;
	acumulador=0;
	resultado=0;
		
		for(int i=0;i<nBarras;i++)
		{
			frecuencia=cvGetReal1D(hist->bins,i);
			acumulador=acumulador+(frecuencia*i);
			nElementos=nElementos+frecuencia;					
			
		}		
		
		resultado = acumulador / nElementos;
		media=resultado;
		nValores=nElementos;
		return resultado;
}

float histograma::obtenerVarianza() //calcular la varianza del histograma
{
	int frecuencia;
	double sumaCuadrados;
	
	sumaCuadrados=0;
	
	if (media==0)
		this->obtenerMedia();
	
	for(int i=0;i<nBarras;i++)
	{
			frecuencia=cvRound(cvGetReal1D(hist->bins,i));
			sumaCuadrados = sumaCuadrados + (frecuencia*(pow((i-media),2)))	;			
			
	}	
	varianza = sumaCuadrados/nValores;
	return varianza;
		
}

float histograma::obtenerDesvEstandar() //calcular la varianza del histograma
{
	if (varianza==0)
		this->obtenerVarianza();
	return sqrtf(varianza);
}

int histograma::obtenerModa() //calcular la moda del histograma
{
	int x;
	float xx;
	
	cvGetMinMaxHistValue( hist, 0,&xx,0, &x);
	return x;
}

float histograma::calcularSesgo(float estimador) //calcular el sesgo del estimador que se pase como parametro
{
	if (media==0)
		this->obtenerMedia();
		
	return (media - estimador);
}

float histograma::obtenerMediana() // calcular la mediana del histograma
{
	vector<int>  valores;
	int frecuencia,n,a,b;
	
	for(int i=0;i<nBarras;i++)
	{
			frecuencia=cvRound(cvGetReal1D(hist->bins,i));
			for(int j=0;j<frecuencia;j++)
			{
				valores.push_back(i);
			}
		
	}
	
	n=valores.size();
	
	if((n%2)!=0)
		return valores.at(((n+1)/2));
	else
	{
		a=valores.at(n/2);
		b=valores.at((n+1)/2);
		return ((a+b)/2);
	}
			
}

int histograma::obtenerRango() //calcular el rango del histograma
{
	
		
	vector<int>  valores;
	int frecuencia,n;
	
	for(int i=0;i<nBarras;i++)
	{
			frecuencia=cvRound(cvGetReal1D(hist->bins,i));
			if (frecuencia!=0)
			{
				for(int j=0;j<frecuencia;j++)
				{
					valores.push_back(i);
				}
			}
		
	}
	n=valores.size();
		
	return (valores.at(n-1)-valores.at(0));
	
}

int histograma::obtenerUmbralOptimo() //metodo OTSU para calcular el valor umbral optimo
{
	if (media==0) this->obtenerMedia();
	
	
	float sum=0,sumB,wB,wF,varMax;
	int frecuencia=0,umbral=0;
	for(int i=0;i<nBarras;i++)
	{
		frecuencia=cvRound(cvGetReal1D(hist->bins,i));
		sum = sum + (i*	frecuencia);	
	}
	
	sumB=0;wB=0;wF=0;varMax=0;frecuencia=0;
	
	for(int i=0;i<nBarras;i++)
	{
		frecuencia=cvRound(cvGetReal1D(hist->bins,i));
		wB=wB+frecuencia;
		if (wB==0) continue;
		wF=nValores-wB;
		if (wF==0) break;
		
		sumB=sumB+(i*frecuencia);
		
		float mB = sumB/wB;
		float mF = (sum-sumB)/wF;
		
		float varBetween=wB*wF*(mB-mF)*(mB-mF);
		if (varBetween>varMax)
		{
			varMax=varBetween;
			umbral=i;
		}
	}
	return umbral;
}




	

